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Abstract 

We present a convergence proof for higher order implementations of the projec¬ 
tive integration method (PI) for a class of deterministic multi-scale systems in 
which fast variables quickly settle on a slow manifold. The error is shown to 
contain contributions associated with the length of the microsolver, the numeri¬ 
cal accuracy of the macrosolver and the distance from the slow manifold caused 
by the combined effect of micro- and macrosolvers, respectively. We also pro¬ 
vide stability conditions for the PI methods under which the fast variables will 
not diverge from the slow manifold. We corroborate our results by numerical 
simulations. 
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1. Introduction 

Many problems in the natural sciences are modelled by multidimensional 
ordinary differential equations with entangled processes running on widely 
separated time scales. One is often interested in resolving the behaviour of 
the slow processes over a long, macro time scale. However, the fast processes 
prevent direct solution of the system by traditional numerical methods. Re¬ 
cently two numerical methods designed to overcome the restriction to the small 
integration time step in these stiff dynamical systems have been much studied; 
the projective integration method within the equation-free framework and the 
heterogeneous multiscale methods (HMM). Each method exists in multiple 
formulations; in the PI method, we mention [1, 2, 3, 4, 5, 6, 7, 8], and in the 
HMM, [9, 10, 11, 12, 13, 14, 15]. There is some debate on the similarities and 
differences between the methods; the interested reader is referred to [16, 17] for 
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a discussion. 

Both methods assume that the fast variables in the full multiscale system 
quickly relax to a slow manifold, after which the dynamics of the slow variables 
is governed by a slow reduced system. Both methods estimate the effective 
influence of the fast variables on the dynamics of the slow variables by 
employing a microsolver to perform short fine-scale computations with small 
time steps (microsteps). This information is used to propagate the dynamics 
on the slow manifold for large time steps (macrosteps) in the macrosolver. 

The philosophy behind each method is slightly different. The PI approach 
estimates the effective slow vector field via direct numerical evaluation, not 
assuming any knowledge on the form of the reduced vector field; this forms 
part of the equation-free approach. In contrast, the HMM philosophy utilises a 
priori analytical knowledge about the reduced vector field. 

In this paper, we focus on numerical methods that are seamless; that is, the 
numerical methods do not explicitly separate the slow variables and the fast 
variables at any stage in the solver, but instead propagate all variables simulta¬ 
neously. These methods are useful in systems where conceptually there exists 
a decomposition or transformation of the system into slow and fast variables, 
but where this transformation is unknown. The added complication of seamless 
numerical methods is that the fast variables are propagated simultaneously 
with the slow variables with the large time step of the macrosolver. This may 
lead to a more severe departure of the fast variables from the slow manifold 
over the macrosteps in comparison to nonseamless methods. 

In first order PI methods the micro- and macrosolver are applied sequentially, 
so the error accrued by the micro- and macrosolver can be analysed separately, 
as for example in [11, 18]. There are two different approaches to extend PI 
to higher order solvers. First, one can still apply the micro- and macrosolver 
sequentially, as in [19, 12, 5, 20]. The analysis in [12, 20] shows that such 
schemes can be accurate to second order in the size of the macrosolver. 
Alternatively, one can apply the microsolver multiple times during each time 
step of the macrosolver, as in [11, 4, 21]. The numerical schemes that we will 
consider take this approach. The analysis of such methods is complicated 
by the requirement that the errors accrued by the micro- and macrosolvers, 
which are intertwined due to the nonlinear nature of the dynamics, have 
to be estimated simultaneously. In [11], an error bound is proposed for a 
seamless HMM scheme of arbitrary order, albeit without proof. In [4, 5, 20], 
second order PI schemes are proposed and analysed. In [21], error bounds for 
the slow variables and stability conditions are derived for an arbitrary order 
Runge-Kutta macrosolver applied to a kinetic equation with linear relaxation. 

In this paper we present a higher order seamless multiscale method as consid¬ 
ered in [11, 4], for a system of nonlinear stiff ordinary differential equations. 
We propose a slight modification of this method which, involving an additional 
application of the microsolver, constructs slow vector fields pointing towards the 
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slow manifold. Both schemes reduce to Runge-Kutta methods if the microsolver 
is switched off. We establish rigorous convergence results for the slow variables 
of these methods. We find that both methods incur error terms propotional to 
the order of the macrosolver, the distance of the fast variables from the slow 
manifold, and an additional term due to the microsolver, independent of the 
order of the microsolver. This result confirms for the two methods we consider 
the error bound suggested in [11]. Furthermore, we find that the error due 
to the microsolver is smaller in our proposed method when both methods are 
employed at the same computational cost. 

A known problem in seamless methods is that the macrosolver may lead 
to a departure of the fast variables from the slow manifold. To combat 
this divergence of the fast variables, several methods have been introduced 
[22, 23, 24, 25]; analytical bounds on the departure of the fast variables from 
the slow manifold over a macrostep have received relatively little attention 
(with the notable exception of [12]). Estimates of the maximal deviation of 
the fast variables from the slow manifold are particularly important when 
bifurcations occur or when the dynamics transits to different solution branches 
(e.g. [19, 1, 7, 26]); if the departure from the slow manifold is too large, the 
transitions may be premature. 

We establish bounds on the departure of the fast variables from the slow 
manifold over the macrosolver. The bounds show that the numerically induced 
departure of the fast variables from the slow manifold scales one order better 
in the macrostep size in our modihed version of PI. Furthermore, these bounds 
allow us to derive stability conditions for both methods under which the 
departure of the fast variables from the slow manifold remains Hnite over the 
macrosteps. 

The paper is organized as follows. In Section 2 we discuss the class of dynam¬ 
ical systems studied, and briefly summarize in Section 3 classic Runge-Kutta 
methods for these systems. We then present two multiscale methods which en¬ 
able the solution of these systems with macro length time steps in Section 4. 
In Section 5, the main part of this work, we derive rigorous error bounds for 
those numerical multiscale methods. In Section 6 we present results from nu¬ 
merical simulations corroborating our analytical findings. We conclude with a 
discussion in Section 7. 

2. Model 

We consider deterministic multiscale systems of the form 

Ze=F{Ze,e) , (2.1) 

with Ze € and time scale separation parameter 0 < £ ^ 1. We assume there 

is a (possibly unknown) decomposition Zi. = {xg,yg) into fast variables XeGM™ 
and slow variables G M" which evolve according to 


ye = g{Xe,ye) , 


( 2 . 2 ) 
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ie=^f{xs,ye)- (2.3) 

£■ 

We consider here the particular fast vector fields of the form 

f{Xe,ye) = ^{-Xe + ho{ye)). (2.4) 

We assume there is a coordinate system such that the matrix is 

diagonal with diagonal entries Xu > 0. We further allow for a scaling of time 
such that min(Aii) = l and define max(Aii) = A. We assume that there exists 
a slow manifold X = hs{y) = ho{y) + O{e), towards which initial conditions are 
attracted exponentially fast. On the slow manifold, the dynamics slows down 
and is approximately determined by 

Y = G{Y), (2.5) 

with Y = ys+0{s) and reduced slow vectorfield 

G{y)=9{he{y),y)- ( 2 . 6 ) 


3. Runge-Kutta Solvers 

We denote by Ze{t"‘) the solution of (2.1) evaluated at the discrete time t” = 
nAt, and by z" the numerical approximation of Ze(t") given by a Runge-Kutta 
solver of order P. Runge-Kutta solvers form approximations to the dynamics in 
terms of increments. For simplicity, we restrict our analysis to Runge-Kutta 
methods in which increments are given recursively by 

kj{z"') = + ajkj-i,£) , (3.1) 

for j = l,2,...,P. The values of the nodes aj depend on the order P (see for 

instance [27]), and satisfy 0<aj<l, with ai=0 so that the first increment is 
defined explicitly. Each increment kj evaluates the vector field P of (2.1) at the 

intermediate time f^ + ajAt. The increments are averaged to define with 

p 

+ (3.2) 


where the weights bj satisfy the condition depend on the order 

P and the particular choice of nodes aj. For instance, for P = 4, the widely used 
fourth-order Runge-Kutta scheme, the nodes and weights may be given by 


b,= 


2 ’ 2 ’ ^ 
1111 
6’ 3’ 3’ 6 


(3.3) 

(3.4) 
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For any P, the nodes and weights are determined such that the application of 
a single Runge-Kutta step of order P to a system with initial condition Ze(t) 
and time step At produces an approximation to ^^(t + At) accurate to within 
see for instance [27]. In particular for linear systems 

A 

Xe = - Xe , 


for which 


d^Xe 

dp 



a single Runge-Kutta step of P-th order can be written as 


X 


n+1 



where the linear amplification factor p is given by the Taylor polynomial to 
order P of an exponential function 


p 


Piv) = Yl 

j=o 


j! 


(3.5) 


A straightforward implementation of Runge-Kutta methods to simulate stiff 
dynamical systems such as (2.2)-(2.4) would be computationally too costly, as 
the time step is restricted to At < 0(e) to ensure numerical stability. 


In the next section we present two numerical multiscale schemes which are de¬ 
signed to overcome the problem of stiffness presented above. These schemes 
employ a microsolver to relax the fast variables towards the slow manifold. 
Utilising the slowness of the dynamics on the slow manifold allows for the ap¬ 
plication of Runge-Kutta methods with large macro time steps At s. 


4. Numerical Multiscale Methods 

We consider two seamless projective integration methods. The first is a general 
order formulation of PI as proposed in [11, 4, 12]. We call this scheme PIl. 
The second is a modification of PIl, which employs information from the mi¬ 
crosolver to define increments which point in the direction of the slow manifold, 
at the cost of one additional application of the microsolver^. We call this 
method PI2. The PIl and PI2 schemes differ in the definition of the increments. 


^We ensure that the overall cost of PIl and PI2 is the same when they are compared 
numerically by adjusting the total number of microsteps in each method (see Section 6). 
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Denote by z" the numerical approximation given by the multiscale scheme to 
Ze{G‘)-, using z” as the initial condition, both methods employ a microsolver 
with small microstep (5t, and then evaluate the vectorfield over a large macrostep 
At^s. Iterating these steps enables one to construct increments which cover 
a macro time scale. The macrosolver then combines these increments in a 
weighted sum in Runge-Kutta fashion. 

We denote by the flow map for the microsolver run for m microsteps with 

time step St and assume that it describes an explicit numerical method of order 
p. We do not specify which particular numerical method is chosen; as we will 
see in Proposition 5.6, increasing the order of the microsolver does not improve 
the predicted overall error scaling. 

In the following we detail PIl and PI2 and highlight their differences. The 
procedures are illustrated in Figure 1 for PII and in Figure 2 for PI2. 

4-1. Projective Integration Scheme PIl 

We describe here a general order formulation of projective integration along 
the lines of [11, 4, 12, 18]. We remark that this formulation is an instance 
where PI and HMM are essentially the same (see [11, 16]). The scheme PIl is 
a modified Runge-Kutta scheme in which the microsolver is employed to relax 
the fast variables close to the slow manifold before each increment is estimated. 


We denote by the approximation of the fast and slow variables at the ru¬ 
th microstep of the j-th increment at time step n, and denote by Mj the integer 
number of microsteps taken before the j-th increment is estimated. We denote 
discrete times associated with microsolvers by subscripts and those associated 
with macrosolvers by superscripts. 

The increments cover a time step of At and are given by evaluating P after an 
application of the microsolver, with 

%(^”) = AtJ^(z](:f^,e) , (4.1) 


for j = 1,2,...,P, where we define for j = 1,2,...,P, m = l,2,...,Mj, as the 
output of the microsolver 



with initial condition 

^ forj = l 

I for j>l 


(4.2) 


(4.3) 


The nodes aj are those used in the increments of a Runge-Kutta solver of order 
P; i.e for P = 4, aj may be given by (3.3). For more general Runge-Kutta solvers 
for PI methods, see [21]. Construction of the microsteps z^^l is illustrated in 
Figures la and Ic, and construction of the increments kj in Figures lb and Id. 
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The macrosolver is then given by the weighted sum 

(4.4) 

j=i 

where the weights bj are appropriate to a Runge-Kutta solver of order P; i.e for 
P = 4, bj may be given by (3.4). The macrosolver is illustrated in Figure le. 
Note that for Mj =0 for all j, i.e. without the microsolver, the scheme reduces 
to a standard Runge-Kutta solver of order P applied to the system (2.1). It is 
not true that PI schemes in general reduce to a numerical discretisation of the 
underlying multi-scale dynamical system if the microsolver is switched off (see 
for example [20]). 

For the analysis of the PIl scheme, it is helpful to explicitly identify the slow and 
fast variables. We therefore decompose the PIl variables into fast and slow 
components (x”,?/"), and into the fast and slow components 

Furthermore, we split the PIl increments kj into fast components and slow 
components kyj, with 

kx,j{x^,y"') 

The macrosolver is then written as 

p 

i=i 

P 

x"^'^^=x'^\+'^bykxj{x”^,y'^) ■ (4.7) 

i=i 

4-2. Projective Integration Scheme PI2 

We present here a modification of the PIl scheme in which the increments 
are given by differences between endpoints of the microsolver. We again denote 
by the approximation of the fast and slow variables at the m-th microstep 
of the j-th increment. The PIl increments are given by (4.1), which we recall 
here as 

= . ( 4 - 8 ) 

for j = 1,2,...,P, where is now defined for j = 1,2,...,P4-1, m = l,2,...,Mj, 
as the output of the microsolver 



= Mg(xli],yP^ , (4.5) 

= ■ ( 4 - 6 ) 


(4.9) 
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with initial condition 




0 



z2i]+ajkj-i{z^) 


for j = 1 
for j > 1 


(4.10) 


Construction of the microsteps is illustrated in Figures 2a, 2c and 2e. The 
PI2 increments are constructed by approximating the vector field T according 
to leading to 

--K) ■ o-ii) 

®i+i ^ ' 

for j = l,2,...,P. The nodes aj with j = l,2,...,F are again those used in the 
increments of a Runge-Kutta solver of order P, and we set op+i = 1. The con¬ 
struction of the PI2 increments kj is illustrated in Figures 2c and 2e. 

Each PI2 increment covers a time step of + We fix the total 

number of microsteps Mj for j > 1 with 


Mj=ajM , 


(4.12) 


for j = 1,2,...,P-|-1 and for some M satisfying ajM SN, so that each increment 
covers a uniform time step of At + M6t=:t^. Note that (4.12) allocates more 
microsteps after larger increments and less after shorter increments. 

The macrosolver is now constructed as a weighted sum over the relaxed incre¬ 
ments kj rather than over kj, with 

p 

(4.13) 

i=i 

where the weights bj again correspond to a Runge-Kutta solver of order P. Note 
that =n{tA +MiSt) for PI2. The macrosolver is illustrated in Figure 2f. 
Again for Mi =M = 0, i.e. without the microsolver, the PI2 scheme reduces to 
a standard Runge-Kutta method of order P. 

As with the Pll scheme, it is helpful to explicitly identify the slow and fast 
variables in the solver. We therefore decompose the PI2 variables z" into 
fast and slow components (a;",?/"), and z^^_ into the fast and slow components 

and we split the PI2 increments kj into fast components kxj and 
slow components kyj, with 

= > (4-14) 

kx,j{x^,y'^) = -^(xlf+^-xlA , (4.15) 

depending via (4.10) on the Pll increments kj. For completeness we recall these 
as 


fcy.a (a;” > y”) = 5 . Vm] ) 


(4.16) 
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= . (4.17) 

The macrosolver is then given by 

i-1 

p 

x"'^^=XM]+'^bjk^j{x^,y'^) . (4.18) 

i=i 

For ease of exposition we write the slow dynamics as 

2/”+'=J/”+ff(x",2/"), 


where the vector field of the slow variables in the PI2 macrosolver, g, is given 

by 

P 

g(a;",2/”) = y)(^|-y" + ^6jfcyj(a;”,y") . (4.19) 

i=i 

Comparing Figures 1 and 2, in the PI2 method the increments point in the 
approximate direction of the slow manifold, so that the macrostep initialises the 
fast variables close to the slow manifold after a macrostep. By comparison, the 
PIl increments can depart from the slow manifold with larger scale separations 
or for initial conditions off the slow manifold. 

5. Error analysis for Projective Integration 

We provide rigorous error bounds for the slow variables of PI in the 
formulations PIl (4.1)-(4.4) and PI2 (4.9)-(4.13), following the general line of 
proof used in [11]. Therein the result for PIl was stated, albeit without explicit 
proof. Furthermore, we establish bounds on the departure of the fast variables 
from the slow manifold over the macrosolver, yielding stability conditions for 
the fast variables. 

Throughout this work we assume the following conditions on the growth and 
smoothness of solutions of our system and on the numerical discretization 
parameters of PI. 

Assumptions 

Al: The zeroth order approximation of the slow manifold ho{y) is Lipschitz 
continuous; that is there exists a constant Lh such that 


\ho{yi)-hoiy2)\<Lh\yi-y2\■ 
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A2: The vectorfield g{x,y) is Lipschitz continuous; that is there exists a con¬ 
stant Lg such that 

\9{xi,yi)-g{x2,y2)\<Lg{\xi-X2\ + \yi-y2\) ■ 

A3: The second order derivatives of ho are all bounded; that is there exists a 
constant Lh' such that 


sup 


l«l=2 


< Lh', 


where we used multi-index notation. 

AJ^: The vectorfield g{x,y) is bounded for all x,y; that is there exists a con¬ 
stant Cg such that 


C's = sup|5(x,2/)|. 

A5: The reduced slow dynamics Y{t) is of class ( 7 max(P,p). jg exist 
constants Cp and C* such that 


C*p = sup 
C* = sup 


d^y(t) 
dtP ’ 
dPY{t) 
dtP ’ 


and in particular there exists a constant CJ satisfying 


Ca =sup|y(t)|. 


A6: The total time At of the macrostep is sufficiently short so that, employing 
the practical constraint M6t<At, 


LQM5t ^ LqAI — . 

Remark 5.1. Assumptions (A1)-(A2) imply that the reduced slow dynamics 
( 2 . 5 ) is also Lipschitz continuous and there exists a constant Lq < Lg(l + Lh) 
such that 


\G{Y,)-G{Y2)\<Lg\Yi-Y2\ . 

Assumption (A4) implies that the reduced slow dynamics is also bounded and 
there exists a constant Cq < Cg such that 


CG = sup|G(y)|. 
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The global Lipschitz conditions can be relaxed to local Lipschitz conditions by 
the usual means. 

We will establish bounds for the error f ” between the PIl and PI2 estimate y" 
and the solution of the full system 


Theorem 5.2 (Convergence). Consider schemes PIl and PI2 run with a 
Runge-Kutta method of order P for the macrosolver and an explicit scheme 
of order p for the microsolver. Given assumptions (Al)-(A6), there exists a 
constant C such that on a fixed time interval T, for each n such that ntA < T, 
the error between the PIl and PI2 estimates and the exact solution of the full 
multiscale system (2.1) are bounded by 

Here a = minj>iaj, p is the linear amplification factor (3.5) for the microsolver 
of order p measuring the attraction of the fast variables to the slow manifold 
over a microstep, and :=max o<i<n, la^o^ — /io(yo ^)| is the maximal de- 

l<k<P+l 

viation of the fast variables from the approximate slow manifold accrued over 
the integration time. 

It is worthwhile to briefly discuss the bound on 5". The term proportional 
to reflects the convergence of the underlying Runge-Kutta numerical scheme 
of order P in the macrosolver. The term proportional to MSt is incurred by the 
drift of the slow variables over the microsteps before estimating the increments 
(regardless of the order p of the microsolver). The terms proportional to the 
time scale parameter e represent the error made by the reduction as well as 
an additional error incurred during the drift of the slow variable over the 
microsteps. The term proportional to (£/<a+ p“'^(—<5t/e))|dmaxl measures 
the mismatch between the slow vector held g{x,y) after an application of the 
microsolver and the reduced vector held G{y). 

We also provide bounds on the deviation |(i"| of the fast variables a;" from the 
slow manifold ho(?/") for PIl and PI2, 

K| = K-ho(2/")|. 


Theorem 5.3 (Stability of the fast variables). Consider schemes PIl and 
PI2 run with a Runge-Kutta method of order P for the macrosolver and a for¬ 
ward Euler scheme for the microsolver. Given assumptions (Al), (A3) and 
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(A4^), the fast variables do not diverge over the macrosolver, so that the largest 
deviation of the fast variables from the slow manifold *■5 finite, if 


\At 

s 



<1 . 


Then the distance of the fast variables from the slow manifold after the n-th 
macrostep satisfies for the PIl scheme the recurrence relation 


i=i 


XAt 

e 


1- 


5t 


aM \ ^ 


\d’^\+LhCg{l + X)At, 


and for the PI2 scheme 


A P 

< ^ Eb, 

a ^ 

j=i 


XAt 


1 - 


St 


iM\ 1 


\d^\+2LyCltl . 


Remark 5.4. The stability condition for the fast variables of the PIl and PI2 
methods is identical to the corresponding stability condition for an Euler macro¬ 
solver given by Assumption 8 in [18] (see also [21]). 


We note that Theorem 5.3 can be formulated for a microsolver of order p>l 
but, as we shall see, optimal convergence results are given by a forward Euler 
microsolver. 


We briefly discuss the stability condition and the bounds for estab¬ 

lished above. The stability condition can be understood as follows: (1 — 6t)e)‘^^ 
denotes the exponential contraction of the fast variables towards the slow 
manifold during the application of the microsolver; if this contraction rate does 
not bring the fast variables within a neighbourhood of e/A of the slow manifold, 
the fast variables will not have sufhciently relaxed and their dynamics remains 
stiff, possibly causing numerical instability over the subsequent integration 
steps. 

The bounds for the deviation of the fast variables from the slow manifold 
are different for PIl and PI2. In particular, Theorem 5.3 suggests that for 
a given macrostep size the fast variables deviate less from the slow mani¬ 
fold in our modified version PI2. This will be confirmed numerically in Section 6. 

In the next section we prove Theorems 5.2 and 5.3. We formulate the proofs for 
PI2 and point out where and how they will differ for PIl. 

5.1. Error Analysis 

We split the error 5" between the PI approximation of the slow variables and 
their true value into two parts. Denote by E(t”) the time-continuous solution 
of the reduced ordinary differential equation (2.5) evaluated at time f”, then 


£- = \ye{n-yl 
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where the first term describes the error between the exact solutions of the full 
system (2.2)-(2.4) and the reduced slow system (2.5), which we label reduction 
error, with 


E^ = y,{e)-y{n, (5.1) 

and the second term the error between PI and the exact solution of the reduced 
slow system (2.5), which we label discretization error, with 

E2 = y--Y{n. (5.2) 

We will bound the two terms separately in the following. 

5.2. Reduction error 

Setting the initial conditions close to the slow manifold with j/£(0) = F(0) + 
and Xg{0) = hg{yg{0)) + we formulate the following theorem for the 
reduction error i?". 

Theorem 5.5. Given assumptions (Al)-(A3), there exists a constant Ci such 
that on a fixed time interval T, for each < T, the difference between the exact 
solutions of the reduced and the full system is bounded by 


\EJf\<Cis, 


with 

C'i=max(|ye(0)-r(0)|,Ly|4(0)|))e^«(i+^'*)‘", 

where de=Xs — he{ye) measures the distance of the fast variables from the slow 
manifold. 

The proof is standard and is omitted here. The interested reader is referred to, 
for example, [28, 18]. 

5.3. Discretization error 

We bound the discretization error E'f =y'^ — Y{t'^) in stages. We hrst give 
a proof for the convergence of a PI approximation of the reduced dynamics 
to the true reduced dynamics in Proposition 5.6. We then compare the PI 
approximations of the reduced and the full multi-scale dynamics, and combine 
the two results to bound Eff. 

To achieve the first bound we introduce the auxiliary vector field G, which 
describes the PI2 method applied to the reduced slow system (2.5). We hrst 
show that GiY (t")) is close to a standard Runge-Kutta solver applied to Y (t” + 
Mi6t); then we bound the difference between the auxiliary vectorheld G{y^) and 
the PI2 vectorheld for the slow variable g(a;”,j/"). 

Denote by the how map for the microsolver of order p applied to the 

reduced system (2.5) for m microsteps with time step 5t. Given initial condition 
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F” at t = t”, we construct G analogously to the construction of g used in PI2. We 
define for j = 1,2,... ,P +1, m= 1,2,... ,Mj as the output of the microsolver 



analogous to (4.9), with initial condition 

yno_[y\ ^ for, = l 


(5.3) 


(5.4) 


analogous to (4.10) and (4.8). The increments are constructed by 

= , ( 5 . 5 ) 

Oj + l ^ / 

analogous to (4.11). Combining (5.3)-(5.5), we form the auxiliary vectorfield 


p 

G{Y-) = Y-f-Y- + Y,b,K,{Y-) , (5.6) 

1=1 

analogous to the PI vectorfield (4.19) of the macrosolver. 


In the following Proposition we demonstrate that G evaluated at Y (t") incurs an 
error of order 0{t^^) over one macrostep, like standard Runge-Kutta methods, 
with an additional error term incurred by the applications of the microsolver. 


Proposition 5.6. Given assumptions (Al), (A2), (A4) and (A5), G{Y{t”‘)) 
provides a numerical estimate of the reduced slow vectorfield with 

Y{r+^)=Y{e)+G{Yir)) + Oit^+\tAM6t), 

where the error term 0{t^^ fiAMdf) is bounded by J-C^tAMSt. 

Proof. The increments Kj of a Runge-Kutta solver of order P applied to the 
reduced system (2.5), initialised at Y with time step tA, are given by 


KfiY)=tAG{Y + a,k,_,{Y)). (5.7) 

For a Runge-Kutta solver of order P initialised at YMi5t) we have 

p 

F(r+i)=F(r + Mi5t)-P^6jK,(F(r + Mi,5t))+0(t£+i) , (5.8) 

1=1 


where the 0{t^^) term is bounded by [27]. Similarly, a microsolver 

of order p satisfies 


r(r+Mi<it)-r”;i 


<c;Mi6tp+^. 


(5.9) 
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The Runge-Kutta solver (5.8) is rewritten as 

y(r+i) +j^b,K,{Y^l) (5.10) 

J-1 

P 

+ (K,{Y{t- + M^5t))-K,{Y^l)) +0 . 

i=i 

Employing assumptions {Al)-{A2) on the Lipschitz continuity of the reduced 
dynamics, (5.9) and the definition (5.7) of the Runge-Kutta increments K we 
bound 

Kj iY{e + Mi6t)) - K, (FjO;')I <LGtA \yP + M^5t) - Y^^ 

+ LatAaj (F(r + M^dt)) - (F^JJ;^) 

<LGtAC;M^5tP+^ 

+ LGtAaj |x,_i(F(r + Mi,5t))-i?,_i(FjO;i) . 
Iterating this relationship with ai = 0 yields 

KjiY{e + MiSt))-Kj(Y;;jH <LGtAC;MiStP+^ +OitlMiStP+^) . 

Upon substitution into (5.10) we obtain 

F(r+i)=F;^f+5]6,X,(F;);i) + (!l(t£+\Mp5tP+i) , (5.11) 

j=i 

which describes a Runge-Kutta method of order P, initialised at Y^’\ The 
auxiliary vectorfield G given by (5.6) with F" = Y(t'^) is now constructed from 
(5.11). We write 

F(r+i) =f;};1 (^K, (F(r)) +K,(F”;i) - K,{Y{e))) +0{tl+\M,5tP+^) 

i-1 

p 

=Yin + GiY{n)+Y,b,{KMf)-K,{Y{n)) +0{t^+\M,6tP+^) , 

i=i 

(5.12) 

where the increments Kj are defined in (5.5). 

We now bound \Kj{Y^’^) — Kj{Y{t'^))\ in (5.12). Rearranging the definition of 
Kj, (5.5), we obtain 
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= + + . (5.13) 

Similarly we use the definition of (5.7), to obtain 

) =^Mi^ +®i+i^AG 

+ aj+iAt +ajKj^^ -GiY^^ YajKj^i)^ 

Y5tG (yY ttjKj—i^ 

= ro”’^+Vaj+iM,5tG(yo”’^'+^) (5.14) 

+ ®i+iAt (g (y^'^ +ajKj_i^ - G(Y^'^ + ajKj_i)^ 

+ a,+^M6t (g (r^;Va,X,_i) - GiY^^’^+^)) , 

where we have suppressed the dependencies of Kj and Kj on the right-hand 
side and used tA = Ait + MSt. Subtracting (5.13) from (5.14), applying absolute 
values and dividing by Oj+i yields the bound 

aj+i \ / 

(5.15) 

+ At g{y^^ YajKj^^ -G{Y^^ +ajKj^i) 

YM5t\G (y”f + a,K,_i) - G(ro"’^+i) . 

We now bound the three lines of (5.15) separately. In the first line, we inter¬ 
pret Yaj+iM5tG{YQ’^^^) as a single Euler step with time step aj+iMSt 

initialised at The remaining term of the first line, , 

describes a microsolver of order p run for Mj+i =aj+iM steps, also initialised 
at Fg”’-^"''^. Therefore 

Fg"’^+Vaj+iM5tG(Fg”’^'+^)-^"^+i’'^‘< 0{aj+iM6tP+\{aj+iM6t)^). 

(5.16) 

The second line in (5.15), employing Assumptions {Al)-{A2) on the Lipschitz 
continuity of the reduced dynamics, is bounded by 

At IG {y^I + - G(F^;Va,Tf,_i) I <a, AtLc \ . (5.17) 

Upon using Assumption {A4) on the boundedness of the reduced dynamics and 
Assumptions (Al)-{A2), we bound the term in absolute values in the third line 
of (5.15) by 

G(F;^f+a,K,_i) -G(Fg"’^+i)| <Lg\y-’^+ 
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—La \ajt\G{Y^^ 

<2G2tA , (5.18) 


where we employed the bound on the nodes aj<l, and Assumption (A5) on 
the smoothness of the reduced dynamics with G 2 =sup|y| =sup|DG(F)G(F)| = 
LgCg- 

Substituting (5.16)-(5.18) into (5.15) yields the bound 




YajlYtLa Kj—i — Kj—i 


+ 2C2tAM5t 


+ 0{M6tP+^,aj+i{MStf)) . 


(5.19) 


Substituting j = 1, noting that oi = 0 and neglecting higher order terms, we have 


Ki(yl)f)-Ki(y(r)) 


<2CMAMSt 


(5.20) 


Iteration of (5.19), seeded with (5.20) at j = l, yields 


K,iY-’^)-K,{Y{e)) 


<20^1 AM5t . 


The Proposition now follows directly by substituting into (5.12) and using the 
weighting condition J2j=i — 

Remark 5.7. Proposition 5.6 can be readily extended for PIl. This Proposi¬ 
tion employs auxiliary increments Kj designed to resemble the PI2 increments 
(4.11). In order to prove this result for the PIl method, one should instead em¬ 
ploy auxiliary increments Kj = AtGfYf^'^), which resemble the PIl increments 

(4.1). Following from (5.12), one can then readily bound \Kj — Kj\ <20^ AtM5t 
to obtain the same bound. 


Proposition 5.6 establishes that using G to propagate the reduced dynamics 
incurs error proportional to t^^ -\-tAM5t. In particular, these terms do not 
depend on the order p of the microsolver. To simplify the calculations, we 
therefore use a forward Euler method as the microsolver for the reduced system. 
In particular, we consider 


m 


= Y-’f,Y5tG{Y-’f,) , 


(5.21) 


for 0 < m < Mj. 

We now use Proposition 5.6 to bound the error between the PI2 approximation 
of the slow variable y” in a full multiscale simulation and the reduced dynamics 

r(G). 
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Lemma 5.8. Given assumptions (A1)-(A5), the discretization error \E'2\ = 
|y” —y(t”)| is bounded by 




pLct" 


g{x\y^)-G{x\y^) 


t\ 


Proof. This result follows from [9]. Employing Proposition 5.6, we have 
E2 =Er 1 - G{Y{e-^)) + 0{tl+\t^M5t) 

=E2-^+Ll-^E2-^+g{x^-\y^-^)-G{y^-^) + 0{tl+\t^M5t) , (5.22) 


where we used the mean value theorem for vector-valued functions to introduce 

£ 2 ;= f DG{Y{r)+e{y^-Y{e)))de, ( 5 . 23 ) 

Jo 

where DG is the Jacobian matrix of G. Recall that the 0(t^'^^,tAMdt) term in 
(5.22) is bounded by Cpt^'^^ YC^tAMSt; taking absolute values of (5.22) then 
yields 


\E-\<{l + \£--^\)\E^ 


n—1 
d 


g{x”-\y”-0-G{y^-0\+C*ptP+^+CltAM6t. 

(5.24) 


To bound |£q ^|, we first obtain an explicit formula for G. Substituting (5.6), 
(5.5) and the Euler microsolver (5.21) into G, (5.6), we obtain 


G{Yn =Y^l - F" + ^ i - F"d) 

j — \ 

Mi-1 


k^O 


M. + i-l 


“j+i V ^ 


Substituting -|-aj_|_iAtG(Y)(^’:^) from (5.4), we obtain the explicit 

formula 

Mi-l ^ h f A^j + i-1 

G(F")=Jt ^ G(F,"’i) + ^^ a,+iAtG(F^)’)+Jt ^ G(F”’^+i) 

fc=0 3 = 1 ^3 + '^ y fc=0 



Substituting into (5.23) with Y^ = Y{t^)+ 6{y^— Y{t'^)), taking absolute val¬ 
ues and employing Assumptions {Al)-{A2) on the Lipshitz constant Lq of the 
reduced dynamics yields 


\CI\< DG(F(r)+%"-F(r))) 

^0 


dO 



John Maclean and Georg A. Gottwald 


19 


Mi-l 


< / |DG(rr') +E-^«.+iAi DG(yi^;/) 

Jo rr; I ^®i+i 


fc=0 

Mt_li — 1 


^ DG(y,">^+') 




d6> 


f(' 


E / I E(5Afx(5t + \ —-— {^LQajj^\lSt-\-LQlVljj^iSi) ] dO . 

J 


Recalling Mj=ajM for j > 1, and using the weighting condition = 

we obtain the bound 


|/:S|<LG(tA+Mi<5t) . 

We substitute this bound into (5.24), with 

+ G|>t£+^ + G2*tAM(5t . 

Iterating the recursive relationship with E^ = 0 yields 

n—1 / 

\E 2 \<J 2 i^ + LG{tA + MiSt)r[c*pt^+^+C^tAM6t 

\gix\y^)-G{y^)\^ 


m—0 


- max 
0 < 2 <n-l 


oLc 


<- 


Lg 


CptX-\-C 2 M 6 t-\- max 
0 < 2 <n-l 


tA 

\g{x\y^)-G{y^)\ ' 

tA 


where =n{tA+MiSt). 

Lemma 5.8 establishes that the error between the PI2 approximation of the 
slow variables of the full multiscale system (2.1) and the true reduced dynamics 
(2.5) contains a term proportional to the order of the macrosolver a term 
proportional to the length of the microsolver M 6 t and an additional term pro¬ 
portional to maxo<i<„_i|g(a;”“^,y”“^) —G(y”“^)|. The latter term measures 
the difference between the PI2 vector field g of the slow variables, and the aux¬ 
iliary vector field G initialised at the same point In order to bound this 

term we define the deviation of the PI approximation of the fast variables from 
the approximate slow manifold over the increments, 

d^^=xZi^-ho{y’;i;^), (5.25) 

for j = l,2,...,P-|-l, m = l,2,...,Mj. The following Lemma bounds 
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Lemma 5.9. Given Assumptions (Al) and (A4), the error between the fast 
variables and the approximate slow manifold during the application of a micro¬ 
solver of order p is bounded for all 0 < m < Mj by 


where 



The first term in the Lemma is a manifestation of the attraction of the fast vari¬ 
ables towards the slow manifold along their stable eigendirection. The second 
term proportional to e describes, as we will see below, the cumulative drift of 
the slow variables y during the microsteps causing a departure from the slow 
manifold for nonconstant ho{y). 

To ensure convergence of the fast variables to the approximate slow manifold 
we require 


0<p 



<1 . 


Proof. Denote the increments of the microsolver as kx,i and ky^i for the fast 
and slow variables respectively, with nodes a' and weights 6'. Employing the 
fast vector field (2.4) we write the fast increments kx,i analogously to (3.1) as 


kxA^:^\y^n=- ^ {x:l=+a'kx,^-l) + +a'rky,^-i) 

= - ^ - ho (C’^) + a'~kx,-i) 


, (5.26) 


where we have used that ky^i=0{5f). Introducing the increment associated 
with a linear system fciin^i((ij)^l) = — AJt -|-a'fciin,i-i) /e, we write 





The linear component /ciiny of the increment can be interpreted as the increment 
of the microsolver applied to the linear system dg = —Adg/e with initial condition 
difh. Therefore, as discussed in Section 3, a microstep taken with the linear 
increments /cuny can be written as 
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Employing (3.2) the microstep is expressed as 


t' 

I / n,j 

‘^m+l ^ lym ) 

i^l 


2 = 1 
P 


=C^+^6'Min,.(C^)+ho(j/;;,’^)+o — 

^ \ p 


2 = 1 


5t^ 


=p(-^)c^+^o(j/r) + o(^) . 


Then 


*m+l 


K-m+l '‘OlS/rri+ll 

ASt 


< 


< 


e 

ASt 




|C^|+L^C,5t + 0( ^ ) . 


+ 0 


St^ 


(5.27) 


The first term in this bound represents the rate of convergence of the fast vari¬ 
ables to the approximate slow manifold; for stability we require |/o(—| < 1- 
The term L^CgSt stems from the drift in the slow variables over a microstep. 
The slowest rate of convergence to the slow manifold is given by min(Aii) = l, 
so we obtain 



<1 . 

Iterating (5.27) then yields 





d 


nj 

m 


LhCgSt + o(^^-^) 


d^^\+Li,CgS + 0{6t) , 


completing the proof of the Lemma. 

Remark 5.10. In the above Lemma we Taylor expand the terms in hgiy^^ + 
a'Jiy^i-i) up to 0{SP/e). However, higher-order terms may improve the error 
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bound. For instance, one can show that for a fourth-order Runge-Kutta miero- 
solver, 


“m+l 


<P 


-j)\d^^ + ^LuC,St 


Remark 5.11. Optimal convergenee of the fast variables to the approximate 
slow manifold during the application of the microsolver is given by a forward 
Euler microsolver with Q<6t< , where 


P 




and the eonvergence rate p^{—5t/e) is bounded above by the exponential con¬ 
vergence exp(—TO(5i/e). The full stability region for the Euler microsolver is 
0 <6t< 2s/X; for further details, see [18]. 


From Lemma 5.9 it follows that the rate of convergence of the fast variables 
to the approximate slow manifold is optimal for an Euler microsolver, and in 
Lemma 5.8 we demonstrated that the dominant error terms between the PI2 
approximation y" and the true reduced dynamics E(t") do not depend on the 
order p of the microsolver. We therefore choose as the microsolver for the PI2 
scheme the forward Euler method to simplify the calculations, and write 


y[d^=y:i^i+stg{x[[f„y[ff,), 


(5.28) 

(5.29) 


We now bound the distance of the PI approximation of the slow variables 
over the microsteps of the full system (2.1), from Yff’G the PI approximation of 
the reduced dynamics over the microsteps. 


Lemma 5.12. Assuming (A1),(A2) and (A6), the PI2 numerical estimate y^if. 
of the slow variable after the application of the microsolver at the j-th increment 
is close to the numerical estimate of the reduced slow variable which was 

initialized at Yff'^ =y'^, with 


('-?)■") 


Ii/m- - I <2Lg 3e + ttj At ( 1 - ^ 


+ e> ( ( Ate + At2(l-^) 


max 

Ji<k<j 


jn,k 


aM 


max 


+ ‘ILqLfiCqCLjt^E 




for 2<j <P, where aM = miuj > i Mj, and 


\y 


n,l 

Ml 




+ 0{MiSte) , 


for j = l. 
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Proof. Employing the definition of the Euler microsolvers (5.28) for the PI2 
scheme and (5.21) for the reduced scheme, and Assumptions {A1)-(A2) on the 
Lipschitz continuity of the reduced dynamics gives 

+ 5t\G{yl,]_,)-G{Y^^^_,) 

<{l + LG5t)\yl,]_,-Y^=_, 

<{l + LG5t)\yl,]_,-Y^^AL,5t\xl,^^^^^ 

=(1 + Lc^f) II + Vt I 

+ Lg6t ho{ylf_-^)-he{y2f^_i) 

={l + LGSt)\y^^^_,-Y;:,’^L,\+Lg6t\dllf^_,\+LgL,e6t + 0{Ste^) , 

where we have defined ho{y) — hs{y) = LgS + 0{s'^). Employing Lemma 5.9 on 
MaZ-iI ^ forward Euler microsolver yields the recursive bound 

te ■ - I < (1 + ) 1_ 1 - Em’/- 11 


+ Lg6t ~ 


d~Lg{LhCgLg)sSt 


which, upon iterating, gives 


<il + LGdtr^ y-’^-Y^‘ 


+ Lg6t Jq’A ^ n-j {l + LG6t)^^~^~ 

k=o ^ ^ ^ 

M,-l 

+ Lg{LiiCg + Lg)sSt {1 + LgSI)^ 


= il + LGStr^ 2/0 ’^-Ec”’'I 


+ L„(5t dn’^ 


(l + LGdt)“^-(l-f)' 


® I ° I Lcdt+f 

, r ^ , r ...A^ + LaSt)^^-! 

Lg[LfiCg-\- L^jsot --——- 

Loot 

LcMjSt 

^ LoMjSt .,n,3 _yn,j , r ^n,j ^ _ 

— ^ yo -^0 ^^9^ “O 
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Lg[LhCg -j- Lg)£dt 


^LcMjSt _ ^ 


Realising that — l<2LaMj5t under Assumption (Ad), we obtain 

\y'lf^-Y;:,’A<2\y^^^-YM+2Lgs\d'S'^\+OieMgSt) . (i 


At j = 1 we initialize at j/q’^ = Yq'^ =2/") obtaining the desired bound 


ylr\-Y^l <2LgS +0{eM,6t) . 


For j > 1, we have 


n,i n,l , 7 

Vo =ym+^jkyj_i 


n.i . A , / n.j—l n,j—l\ 

=2/Mi ) , 


using the definitions (4.10) and (4.16), and 


Y-^^=Y-^^^+agAtG{Y-l-^), 

using (5.4). Substituting these into (5.30) and employing assumptions (Al)- 
(A2), we obtain 


<2 ylil-YM^ +2a,Ai -G(r-r/) 


"d-l , "J'-I'l 


+ 2Lg£: +OiMSte) 


<4Lge|d”’ ^+2ajAt^g{xlil_^ J - J | 

+ 2a, At I Giy ^^-^) - GiY^^'^^^) | + 2LgS |d”’^' | + 0(Mdte,Midfe) 
<JLge Id^’i I + 2LgagAt Id);^-^ I + 2LGagAt\y-^-^ - 
+ 2Lg£ dg’-’ +0{Mdt£,MiSt£) 

I I / dt \ ■ I 

<4Lg£ dg’^ +2LgajAt I 1-I dg’-’”^ +2LgL;iGgajAt£ 


+ 2LGagAt yli^Zl-Y;^^-^ +2Lg£ d”’^ +0(Mdte,Midte) 


where we have employed (5.31) to bound y'Zl — Y^’^ and Lemma 5.9 to bound 
Rearranging and taking the maximum over all increments in the terms 
in Idg’-’l and |dg’^| produces 


max dg’ +2LgL;jGgaAA£ 


Vm]-Ym- ^ 2Lg 3£ + a, At 1 - 
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n,j-l _Yn,j-l 
yM,_i ^M,_1 


-0{MSts,Mi6te) . 


Iterating this relation yields to lowest order 

aM\ 


yMi ^ Mi 


<2lJ 3£ + a, Atf 1- — ] 


max 




+ O Ate + (1 - ^^ 


max 


7n,fc 


+ 2LgLfiCg(Xjt^S 

t\e,M6te,Mi6te ) . 


Lemma 5.12 provides bounds on the difference between solutions of the PI ap¬ 
proximation of the slow variables in the full multiscale system and those of the 
PI approximation of the reduced system during the application of the micro¬ 
solver. We use this result to bound the difference between the vectorfield of the 
PI2 method g given by (4.19) and the auxiliary vectorfield G given by (5.6). 


Lemma 5.13. Assuming (Al)-(A6), the auxiliary vectorfield G is close to the 
vectorfield g with 


\gix-,yn-Giy-)\<2Lj-e + At[ 1 -- 


aM\ 


max \dQ'’^\+2L„LhCntA£ 
l<k<P+l ^ a a 


-0{Mi6te) 


Proof. We write 


ISiP.y-') -G(r)l = - r”;‘ + '£;^ UCl “< “ 

j=l “j + 1 ^ ^ 


<( --1 
a 


n,l A^n.l 
VMi ~ ^Mi 


P , 


, nj + 1 _ + i 

yMj+i ^j+i 


where we have used that 1 < l/a^ < 1/a for j > 1. Employing Lemma 5.12 yields 


|5(x",y")-G(2/”)|<(--l)f2L,£ 


7",1 


+ 2LqLhC„Mi5te 


i=i \ 


Uj+i 

2LgL}iCgtA£ 


A / , dt 

+ At[ 1- 

e 


6tV 


V 


M^ 


max 
f l<k<j 


<2L„ [ (-— 1)£-|-At f 1 — —^ I max Idn’^l 

\ a \ e J I i<k<p+i 

-\-2LgLfiCgtA£~h2(\ - )LgLfiCgM\dt£ . 
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Remark 5.14. In order to prove the above result for the PIl method, one con¬ 
structs G with the increments given by vector field evaluations as discussed in 
Remark 5.7. The Lemma then follows along the same lines as the proofs in 
[11, 18]. 

We are now in the position to establish the bound on the discretization error 

\E^\ = \y--Y{n\, 

which we formulate in the following theorem. 

Theorem 5.15. Given assumptions (Al)-(A6), there exists a constant C such 
that on a fixed time interval T, for each nAt < T, the error between the solution 
of the projective integration scheme PI2 and the exact solutions of the reduced 
system is bounded by 

\E[f\<C M6t-\- e ' ^Mmaxl+£^ > 

where Idmaxl :=max o<i<n, the maximal deviation of the fast variables 

l<k<P+l 

from the approximate slow manifold over the increments and macrosteps. 


Proof. Combining the bound on the discretization error obtained in Proposi¬ 
tion 5.6, 


r>Lrit^ 


\E2\<- 


Lq 


C*pt\+C[M5t+ max 

0<i<n-l 


g{x\y^)-G{x\y^) 

tA 


with Lemma 5.13 we obtain 




\E2\<- 


Lg 


C*pt^ + C^M6t + 2Lg{- — +e 
a tA 


MmaxI 


Theorem 5.2 now follows from Theorems 5.5 and 5.15. 


Remark 5.16. Following the comments in Remarks 5.7 and 5.14, one can 
obtain the same bound for the PIl method as obtained in Theorem 5.15 for the 
PI2 method. 


Besides the parameters used in the numerical scheme, i.e. the macrostep size 
tA, the number of microsteps M with microstep size 5t, and the time scale 
parameter e, the error bound also involves the maximal deviation of the fast 
variables from the approximate slow manifold |(imaxl- 
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We now establish Theorem 5.3 by bounding the distance of the fast variables 
from the slow manifold over the increments and macrosteps in the PIl and PI2 
schemes. This provides stability conditions for the fast variables in the PIl 
and PI2 schemes, i.e. conditions under which i® finite. These stability 

conditions are crucial for the successful application of the seamless PI methods, 
since - as we shall see - the fast variables depart from the slow manifold at rate 
proportional to AAt/e^l over the increments and macrosteps. 

5.4- Stability of the fast variables 

We first bound the distance of the fast variables from the slow manifold over 
the PIl increments kj, which are employed in both PIl and PI2. 


Lemma 5.17. Given assumptions (Al) and (A4), the distance of the fast vari¬ 
ables from the approximate slow manifold after the j-th increment of the n-th 
macrostep in the PIl and PI2 methods, given by = —h(i{yQ’^)\, satis¬ 

fies 



where =(i" =a;" —and where we define aM = minj Mj for the PIl 
method to preserve the notation for both methods. 


The first term in this result measures the combined effect of the convergence of 
the fast variables towards the slow manifold over the microsteps, proportional 
to {l — St/e)°'^, and the departure of the fast variables from the slow manifold 
over the subsequent increment, proportional to XAt/e. 

Proof. The PIl and PI2 methods (cf (4.3) and (4.10)), respectively satisfy 

+aj+ifca;.j(a;”,y") 

=x'fiil+aj+iAtj{-x'lf^+hoiyMl)) ■ 

Then 


jnj + l 
*0 


-hoivo 


A 


+ aj+iAt ^ii-xlf^+ho[y'lf]))-ho (%’^"^^) 


-aj+iAt- (Td + /lo (dMi) - ho (% 


n,J + l' 


<- 


XAt 


wo 

“m,- 


+ Lh 


n,l 71,7 + 1 

Vm^-Vo 


^Mi 


(5.32) 
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Employing Lemma 5.9 with the Euler microsolver and substituting the slow 
component of (4.3) and the definition of k, (4.5) and (4.16) respectively, for 
we obtain 

\ / St \ 

LfiCgS \ LhCgAt-\- il —j \d^\-\- LfiCgS 

|C'|+L/,Cg(l + A)At+h-^j \d^\+LhCge. 
Iterating this recursive relation concludes the Lemma. 

We now prove Theorem 5.3 in two parts. We first establish bounds on the 
distance of the fast variables from the slow manifold after one macrostep of the 
PIl scheme in Lemma 5.18, and then follow with the corresponding bound for 
the PI2 scheme in Lemma 5.20. 



Lemma 5.18. Given assumptions (Al) and (A4), the distance of the fast vari¬ 
ables from the approximate slow manifold after the n-th maerostep in the PIl 
seheme, given by —/io(2/”~''^)|, satisfies the recurrence relation 



In particular, the fast variables do not diverge if 


XAt 

e 



When this condition is satisfied, the distance of the fast variables from the slow 
manifold after a macrostep can be written to lowest order as 

) \d-\ + uCg{i+x)At. 

Proof. We substitute the macrosolver (4.7) and increments (4.6) into 
obtaining 

|d«+i| = |x"+i-/ro(j/"+i)| 

A At ^ 

I n,l \ ^ u u ( n+l\| 

— \^Mi ~ / ~ ^o(y )| 

i=i 
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\/\f ^ 

= I H dn] + ho {VmI ) - ^0 (y”+ ^) I 

i=i 

p 

i=i 

Employing Lemma 5.9, Assumption (Al) on the Lipschitz continuity of the 
approximate slow manifold and (4.3)-(4.5) to bound \ho{y]^^J — ho{y"''^^)\ pro¬ 
duces 

I <^ E 6, (l - f I ^j 

/ ^f\ 

d~Lh yMi~y^^^ \ d^hCgs 

K^\+LnCgil + X)At 
H-fl——J \d'^\-\-LfiCgE . 

Substituting the bound from Lemma 5.17 into Idp’^l yields 





completing the Lemma. 

Remark 5.19. The bound presented above for PI 1 is not sharp. If the duration 
of the microsolver is sufficiently large with M6t^e, or if the fast variables are 
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initialised on the slow manifold with |(i° = 0|, then the PIl increments may be 
approximately tangent to the slow manifold and higher order accuracy in At 
can be achieved. 


We now formulate analogous results for PI2. 

Lemma 5.20. Given assumptions (Al), (AS) and (A4), the distance of the fast 
variables from the approximate slow manifold after the n-th macrostep in the 
PI2 scheme, given by — ho{y'^~^^)\, satisfies the recurrence relation 


(l-—) 

jn+l I ^ V e / 








+ LhCg{l + \)AtY, I 

i -0 V 

A/ \ 

+ 0\e,[l-^j\ \d 


J id’ 

aM\ ^ 


+ 2 Lf,lCgtj^ 


In particular, the fast variables do not diverge if 


V " 


aM\ ^ 




e 


When this condition is satisfied, the distance of the fast variables from the slow 
manifold after a macrostep can be written to lowest order as 




P /xA, / c.\aM\t 




i=i \ 


/ XAt f St 


\ e 


\d^\+2Lh'Chl . 


We remark that the first bound presented in the above Lemma for PI2 is pre¬ 
cisely (l— y) /a times the bound presented for PIl in Lemma 5.18, with 
an additional term proportional to and the curvature of the slow manifold, 
measured by Lh'. 

Proof. We reformulate employing (4.18) and (4.15) and Lemma 5.9 to 
estimate 


|d"+i| = |x"+'-ho(|/"+i) 
P , 

, n.l . 

= X 


mEE ^ - ^mI) - ^o(y”+^) I 

j — 1 

p 


< 


^0 (j/E) “ ^0 (l/”+^) + E ( ^0 ) - do {yfdl )) 

j=i ^ ^ 
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p 




a^+i 


d 


nj + l 

Mj + i 




(5.34) 


The first term on the right-hand side of (5.34) can be Taylor expanded to second 
order to obtain 




i=i 


=-m< 


( E 


bjky 




- E (jlb.ky,) +0{tl) , (5.35) 

\a\=2 \j=l / 


where we used multi-index notation to denote the second order derivatives of 
hp. Similarly, the second term on the right-hand side of (5.34) can be estimated 
by Taylor expanding the chord ho( 2 /)()'’^) —/io( 2 /)()^) to second order, employing 
(4.14), with 




—aj+iD/io(j/^i )fcj/j- 

E + ^(^a) ■ (5.36) 

| a |=2 


Substituting (5.36) and (5.35) into (5.34) yields 

p 




| a |=2 


max 

| a |=2 


E +1 .i I E d 


i=i 


w=i 




i=i 


aj+i 


<2Lh' max |fcyj(x”,y”)|^ 

|q|=2. 

l<fc<P 


+ EE“ (MM^ti^l + MMil) +^(^a) 
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where we used that employed Assumption {A3) on the 

Lipshitz continuity of the Jacobian Dho- Employing Lemma 5.9 and recalling 
that the time step covered by each PI2 increment is Ia, we obtain 






Qj + i 


1-^ 

£ 


aM 




which on substituting Lemma 5.17 becomes 


|d"+i|< 


(1-f) 


aM 


. 1=1 

+ (i- 


\At 


1 - 


St 


aM\ 3 


|d"|+ LhCg{l + X)At 




aM\ n 


-\-2LhiCgt^-\-0{e) . 

The upper bound on |(i"| diverges as n increases unless 


fl_« 


aM 


XAt 

e 


1-^ 

£ 


iM\ ^ 


<1 Vj, 


completing the Lemma. 

Combining Lemmas 5.18 and 5.20 yields Theorem 5.3. 


6. Numerics 


We now illustrate the key results of Theorem 5.2 with a fourth-order Runge- 
Kutta macrosolver, which we recall here for P = A including the constants ob¬ 
tained in the proof. We employ a forward Euler microsolver with p=l unless 
otherwise stated. 

In order to compare PIl and P12 at the same computational cost, we choose 
the number of microsteps in the P12 method proportionally lower so that the 
two methods take the same number of microsteps over one macrostep, with 
Mj =M in PIl and Mj ={M,M/2,M/2,M,M} in PI2. Recalling Theorem 5.2 
for P = 4, the discretisation error |i?^| = |j/" — F(t”)| is bounded in the PIl and 
PI2 method by 






C^t\-\-C2Mj iidt-\-2Lg ( 9-—he 


K 


\+2LgL 


hCg£ 


with M[ = M for PIl and Mjj = M/2 for PI2, and where for P = 4, a = 1/2. 
However, the distance of the fast variables from the slow manifold after a 
macrostep scales differently in the two methods. Recalling Theorem 5.3, 
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the distance of the fast variables from the approximate slow manifold |c?"| = 
|a;” —ft-o(2/”)| is bounded for stable applications of the PIl method by 


\d 


rn+l I 


/AAt/ 5t 


MN 


\d^\+LhCg{\ + \)^t 


and for stable applications of the PI2 method by 




< 2 




\d'^\+2Lh'Clti . 


We show results for the multiscale system 

ye = -Xeye-Oiyl 
. _ -a;£ + sin^(7/e) 


( 6 . 1 ) 

( 6 . 2 ) 


which has stable fixed point at (0,0). At lowest order in £, the associated slow 
reduced system is given by 


Y = -Ysiii^{Y)-aY‘^. 


(6.3) 


For higher order approximations of the slow manifold and the associated 
coordinate transformations relating y and Y the reader is referred to the useful 
webtool [29] (see also [30]). 

The system (6.1)-(6.2) with initial conditions y£(0)>0 is locally Lips- 
chitz with Lipschitz constant Lfi = l and Lg =max(|a;£|+2a|j/£|) where 
the maximum is taken over the local region around the initial condi¬ 
tions (a:£(0),ye(0)) under consideration. The vectorfield of the slow dy¬ 
namics (6.1) is locally bounded by Cg = max(|a:eye|-|-a|yep), with the 
maximum taken over the same region. Note that the free parameter a 
controls the constants CJ = Q;ye(0)^(2a-|-sin(2j/e(0)))-|-(!l(sin^(ye(0))) and 
Q = 16a3ye(0)3-b8a'‘ye(0)®-fC>(sin®(j/£(0))). 


We first investigate how the discretization error \E2\ scales with the macrostep 
size At in the PIl and PI2 methods, when all other parameters are kept fixed 
(except n, to fix the final time T). Our analytical result predicts that, so long 
as is small and the practical assumption e<MSt is satisfied, results 

will be divided into two regimes: for Cpt^KC^MSt, the bound for \EJ^\ is 
dominated by the term proportional to C^MSt and \E^\ is independant of 
for Cpt^> C^MSt, the scaling is \E2\^t^. The slight advantage of the PI2 
method in this case is that distributing the same total number of microsteps 
over more applications of the microsolver results in lower error due to Mdt. To 
keep the term proportional to Idmaxl small in both cases, we choose parameters 
so that (l — < 1- The predicted regimes are clearly visible in Figure 

3, where results are presented for a range of macrostep sizes At for PIl and 
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PI2. We choose a = 0.2, and the scale separation parameter e=10“®. We use 
M = 40 microsteps with microstep size St = OAs, while the number of iterations 
n vary from 20 to 10® to keep T = 1 fixed for all values of At. Initial conditions 
are chosen to lie on the approximate slow manifold with y® = l, =sin^(l). 
The Lipschitz constants are Lg = l.l and L/j = l, the bound on the vector field 
of the slow dynamics is Cg = 2, and the maximal derivatives of the reduced 
slow dynamics are Cl =4 and C| =8. 


We present results for the error scaling of \E2\ with the microstep size 6t 
in Figure 4. To focus on the scaling with M6t, we select parameters with 
Cpt^KC^Mdt, and control the distance of the fast variables from the slow 
manifold Id^axl by ensuring (l — < 1. Figure 4 confirms our analytical 

result, that under the condition Cpt£<C|Mi5t, the discretization error scales 
like E^^MSt. The advantage of the PI2 method here is that distributing the 
microsteps over an additional application of the microsolver leads to an overall 
smaller error compared to PIl due to the smaller drift of the slow variables 
over the microsolver. 

We use a second-order Runge-Kutta microsolver (i.e. p = 2), to demonstrate 
that the scaling with 6t is not affected by the order of the microsolver. We 
choose a=l, and the scale separation parameter £=10“®. We use M = 100 
microsteps and n = 50 iterations of each method. The macrostep size At varies 
from 0.0035 to 0.0032 to keep fixed as 6t increases. Initial conditions are 
chosen to lie on the approximate slow manifold with y® = 5, x** = sin^(5). The 
Lipschitz constants are Lg = ll and Lh = l, the bound on the vector field of 
the slow dynamics is Cg =30, and the maximal derivatives of the reduced slow 
dynamics are C| = 50 and Cl = 2000. 


We illustrate the linear scaling of \Ell\ with the maximal distance jd^axl of 
the fast variable from the approximate slow manifold after a macrostep in 
Figure 5. We do so by scaling the initial condition for the fast variables, x®. 
To ensure that the error is not dominated by the initial initialization error 
we choose parameters which render the scheme unstable, allowing for 
divergence of the fast variables from the slow manifold over the macrosteps, 
i.e. |fi”|>|d°|. Figure 5 confirms clearly the linear dependence of \E'l\ on 
MmaxI- We choose a = l, and the scale separation parameter £ = 10“^. We 
use M = 100 microsteps with microstep size (lt = 0.0l£, and n = 5 iterations 
of each method with macrostep size At =10“®. Initial conditions are ip = 1, 
X® e [sin^(l)-|-0.01,sin^(l)-I-1]. The Lipschitz constants are Lg=3 and Lh = l, 
the bound on the vector field of the slow dynamics is Cg = 2, and the maximal 
derivatives of the reduced slow dynamics are Cl = 3 and Cl = 24. 


We investigate how fbe maximal deviation of the fast variables from 

the slow manifold, scales with At in the PII and PI2 methods. We choose 
parameters satisfying <Cl, so that the bounds presented for \d^\ 
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in Theorem 5.3 imply 

ICaxl<i/.C7,(l + A)At + 0(At2) 

for the PIl method, and 

|dS.axl ^^ At) 

for the PI2 method. As noted in Remark 5.19, the bound for the PIl method 
is not tight for systems with |fi°| =0. We therefore choose initial conditions off 
the slow manifold. Furthermore, to ensure that the initial error \d^\ does not 
dominate the error Id^axl) record Idmaxl after the first macrostep. Figure 6 
clearly shows the linear dependence of |dmaxl with the macrostep size At for 
PIl, and the quadratic dependence of Id^axl with the macrostep Ia for PI2. 
We choose again 0 = 0.2, and the scale separation parameter e=10“®. We use 
M = 40 microsteps with microstep size 6t = 0.4:S, while the number of iterations 
n vary from 20 to 10® to keep T = 1 fixed for all values of At. Initial conditions 
are y^ = l, =sin^(l) +1. The Lipschitz constants are Lg = 5 and Lh = l, the 
bound on the vector field of the slow dynamics is Cg = 2, and the maximal 
derivatives of the reduced slow dynamics are Cl =4 and C| = 8. 

Finally, we investigate how Aj/C'^* = |y7’,At _yT,At/ 2 | the 

macrostep size At where the outputs of the PIl or PI2 methods 

with macrostep size At and final time T. In [4] Ay^’^* was used as a measure 
of the numerical error. In Figure 7 we show how Ay^’^* scales with At 
for a fourth-order Runge-Kutta macrosolver, employing the same numerical 
parameters as in Figure 3. It is seen that — y^At/2| ^ fQj. values of 

At whereas the actual discretization error is dominated by MSt at the lower 
values of At (cf. Figure 3). Hence, such proxies for the numerical error have to 
be treated with caution when evaluating PI methods. 

We comment that the numerical results presented here are robust; in par¬ 
ticular, we confirm that identical scalings can be produced from simulations of 
the Michaelis-Menten system employed in [22], although its fast dynamics does 
not satisfy our form (2.4), and the Brusselator with rapidly replenished source 
employed in [2], where the approximate slow manifold is constant. 

7. Discussion 

We have introduced PI2, a seamless numerical multiscale method with 
a higher-order macrosolver, which is a slight modification of a standard 
implementation of a projective integration method, PIl, involving an additional 
application of the microsolver. In both PIl and PI2, each increment is rooted 
on the slow manifold. In PIl the increments typically do not end on the slow 
manifold. In contrast, the additional application of the microsolver assures that 
in PI2 each increment also ends on the slow manifold, even for slow manifolds 
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with non-vanishing curvature (see Figures 1 and 2). If the slow manifold is 
sufficiently linear over the course of one macrostep, the increments of PI2 then 
all lie approximately tangential to it. 

We presented error bounds for the slow variables for both methods, expressed 
in Theorem 5.2. The error bounds are not affected by the order of the micro¬ 
solver used (though strictly speaking, we only considered explicit microsolver 
schemes). Hence the contribution of the microsolver to the error constitutes a 
bottleneck for PI methods, after which the error in the slow variables cannot 
be improved by adjusting the macrostep size or the order of the macro- or 
microsolver. Hence there is no gain to be expected in the slow dynamics when 
microsolvers other than forward Euler schemes are used. 

In Theorem 5.3 we derived bounds for the unphysical deviation of the fast 
variables from the slow manifold, which may cause numerical instability [12], 
and provided a stability criterion for the macrostep size. 

The Theorems now allow us to compare the PIl and PI2 schemes. A fair com¬ 
parison requires that both schemes are operated at the same computational 
cost. Hence, PI2 utilises less microsteps per application of the microsolver dur¬ 
ing the construction of the increments as the total number of microsteps is 
distributed over one more application of the microsolver. Consequently, the ab¬ 
solute discretisation error of PI2 is smaller when compared to PIl. Theorem 5.3 
establishes that the PI2 method incurs less deviation from the slow manifold 
as the deviations scale quadratically with the macrostep size rather than lin¬ 
early as for PIl. The improved stability can be attributed to the increments of 
PI2 pointing towards the slow manifold, enforced by the additional relaxation 
towards the slow manifold when constructing the increments. 
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(c) (d) 



Figure 1: Sketch of the PIl scheme for a second order Runge-Kutta macrosolver. 
The microsolver is employed in la, Ic, the increments kj are given by vectorfield 
evaluations in lb, Id, and the macrosolver is illustrated in le. For this scheme 
oi = 0, 02 = 1 and 5i = &2 = 1/2. 
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(C) (d) 




Figure 2: Sketch of the PI2 scheme for a second order Runge-Kutta macrosolver. 
The microsolver is employed in 2a, 2c, 2e and the (now auxiliary) quantities kj 
are given by vector field evaluations in 2b, 2d. The increments kj are estimated 
in 2c, 2e, and the macrosolver is illustrated in 2f. For this scheme ai = 0, 
02 = 03 = 1 and bi =62 = 1 / 2 . 
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Figure 3: Plot of \og\E'^\ versus log(At) for fixed time of integration T = 1 of 
the system (6.1)-(6.2). The crosses represent results from the PIl scheme and 
the circles represent results from the PI2 scheme. The dashed line is a linear 
regression line with a slope of 3.93. 
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Figure 4: Plot of logjif^l versus \og{6t) for fixed time of integration T = 0.18 of 
the system (6.1)-(6.2). The crosses represent results from the PIl scheme and 
the circles represent results from the PI2 scheme. The dashed lines are linear 
regression lines with a slopes of 0.99 and 1.0, respectively. 
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Figure 5: Plot of loglFl^'l versus log|d°| for fixed time of integration T = 0.0056 
of the system (6.1)-(6.2). The crosses represent results from the PIl scheme and 
the circles represent results from the PI2 scheme. The dashed lines are linear 
regression lines with slopes of 1.01. 
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Figure 6: Plot of log|c?”a^xl versus log(At) for fixed time of integration T = 1 of 
the system (6.1)-(6.2). The crosses represent results from the PIl scheme and 
the circles represent results from the PI2 scheme. The dashed lines are linear 
regression lines with a slopes of 0.98 and 1.94, respectively. 
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Figure 7: Plot of log|A?/^’^‘| versus log|At| for fixed time of integration r = 1 
of the system (6.1)-(6.2). The crosses represent results from the PIl scheme and 
the circles represent results from the PI2 scheme. The dashed line is a linear 
regression line (of the PIl output) with slope of 4.05. 




